Spis treści

  1. Streszczenie analizy
  2. Wczytanie i przygotowanie danych
  3. Podstawowe informacje o zbiorze danych
  4. Analiza wartości atrybutów
    1. Długość złowionego śledzia (length)
    2. Dostępność planktonu (cfin)
    3. Dostępność planktonu (chel)
    4. Dostępność planktonu (lcop)
    5. Natężenie połowów w regionie (fbar)
    6. Roczny narybek (recr)
    7. Łączne roczne natężenie połowów w regionie (cumf)
    8. Łączna liczba ryb złowionych w ramach połowu (totaln)
    9. Temperatura przy powierzchni wody (sst)
    10. Poziom zasolenia wody (sal)
    11. Miesiąc połowu (xmonth)
    12. Oscylacja północnoatlantycka (nao)
  5. Korelacje między zmiennymi
  6. Zmiana rozmiaru śledzi w czasie
  7. Regresor przewidujący rozmiar śledzia
    1. Przygotowanie zbioru trenującego i testowego
    2. Regresja liniowa
    3. Regresja liniowa z regularyzacją (Ridge, Lasso i Elastic-Net)
    4. Regresja nieliniowa (wielomianowa)
    5. [Gradient Boosting Machine] (#gbm_model)
  8. Analiza ważności atrybutów najlepszego modelu (Gradient Boosting Machine)

Streszczenie analizy

Celem projektu było określenie jakie mogą być główne przyczyny stopniowego zmniejszania się długości śledzi oceanicznych wyławianych w Europie. Aby odpowiedzieć na to pytanie w naszym raporcie dokonałyśmy szeregu analiz mających nas przybliżyć do rozwiązania. Zaczęłyśmy od wczytania i przygotowania danych w którym zadbałyśmy o poprawne obsłużenie brakujących wartości. Następnie dokonałyśmy wstępnego przeglądu danych aby dowiedzieć się podstawowych informacji o rozkładzie wartości poszczególnych atrybutów. W kolejnej sekcji przeanalizowałyśmy możliwe zależności występujące pomiędzy zmiennymi w zbiorze danych. W tym celu zbudowałyśmy macierz korelacji każdej zmiennej z każdą oraz dla wybranych, najbardziej istotnych zależności przedstawiłyśmy je na wykresie w celu zaprezentowania ich charakterystyki. Następnie przedstawiłyśmy interaktywny wykres z użyciem biblioteki plotly, przedstawiający jak zmieniała się długość śledzia w czasie. Na koniec podjęłyśmy się zbudowania regresora przewidującego długość śledzia w czasie. Aby uzyskać jak najlepszą jakość predykcji przetestowałyśmy kilka modeli regresji. Najpierw jednak podzieliłyśmy zbiór danych na dane treningowe, walidacyjne i testowe aby uzyskiwane przez modele wyniki były miarodajne. Predykcje każdego z modeli zostały ocenione za pomocą miary RMSE oraz współczynnika determinacji R-squared.

Zaczęłyśmy od najbardziej podstawowej regresji liniowej. Zgodnie z oczekiwaniami model ten nie był w stanie dopasować się w dobrym stopniu. W związku z tym zdecydowałyśmy się spróbować użyć różnych technik regularyzacji. Co prawda poprawiło to wynik ale tylko nieznacznie. Dlatego wyciągnęłyśmy wniosek że zależności które występują w zbiorze danych są bardziej niż liniowo skomplikowane. Aby uzyskać model lepiej dopasowujący się do danych zbudowałyśmy model regresji wielomianowej poprzez proste dodanie istniejących zmiennych podniesionych do potęgi do modelu. Zabieg ten spowodował znaczny spadek błędu średniokwadratowego a zatem lepsze dopasowanie. Jednak ze względu na prostotę tego rozwiązania postanowiłśmy spróbować skorzystać z algorytmu dedykowanego do tego tego typu problemów, a mianowicie z GBM. Okazał się on być najlepszy z pośród przez jas przetestowanych i zapewnił nam wynik RMSE równy … oraz R-squared równe ….

Wczytanie i przygotowanie danych 

##  [1] "lime"      "gbm"       "glmnet"    "Matrix"    "plotly"   
##  [6] "ggplot2"   "dplyr"     "corrplot"  "stats"     "graphics" 
## [11] "grDevices" "utils"     "datasets"  "methods"   "base"

Kod zapewniający powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych:

set.seed(13)

Wczytanie danych i podgląd postawowych informacji na temat zbioru danych:

raw_data <- read.csv("sledzie.csv")
raw_data <- raw_data %>% 
              rename(t = X)

all_rows_number <- nrow(raw_data)
str(raw_data)
## 'data.frame':    52582 obs. of  16 variables:
##  $ t     : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ length: num  23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
##  $ cfin1 : Factor w/ 40 levels "?","0","0.01",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ cfin2 : Factor w/ 49 levels "?","0","0.01",..: 14 14 14 14 14 14 14 14 14 14 ...
##  $ chel1 : Factor w/ 49 levels "?","0","0.2287",..: 20 20 20 20 20 20 20 20 20 20 ...
##  $ chel2 : Factor w/ 52 levels "?","10.10963",..: 1 23 23 23 23 23 23 23 23 23 ...
##  $ lcop1 : Factor w/ 49 levels "?","0.30741",..: 20 20 20 20 20 20 20 20 20 20 ...
##  $ lcop2 : Factor w/ 52 levels "?","10.72889",..: 23 23 23 23 23 1 23 23 23 23 ...
##  $ fbar  : num  0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
##  $ recr  : int  482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
##  $ cumf  : num  0.306 0.306 0.306 0.306 0.306 ...
##  $ totaln: num  267381 267381 267381 267381 267381 ...
##  $ sst   : Factor w/ 52 levels "?","12.7690663857",..: 38 38 38 38 38 38 38 38 38 38 ...
##  $ sal   : num  35.5 35.5 35.5 35.5 35.5 ...
##  $ xmonth: int  7 7 7 7 7 7 7 7 7 7 ...
##  $ nao   : num  2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...

Z rezultatu wywołania funkcji str wynika, że niektóre kolumny zawierają znaki zapytania, co oznacza brakujące dane.

raw_data[raw_data == "?"] <- NA
missing_rows <- raw_data[!complete.cases(raw_data),]
data <- raw_data[complete.cases(raw_data),]

Przetwarzanie brakujących wartości:

Każdą brakującą wartość, oryginalnie reprezentowaną jako znak zapytania, zamieniamy na charakterystyczne dla języka R wyrażenie NA. Dzięki temu możliwe jest wykorzystanie funkcji complete_cases.

Podsumowanie liczby wierszy w zbiorze danych przed i po pominięciu wierszy z brakującymi danymi:

## [1] "Przed: 52582, po usunięciu wierszy zawierających '?': 42488"

Podstawowe informacje o zbiorze danych

Zbiór danych powstał na potrzeby analizy rozmiaru śledzi oceanicznich wyławianych w Europie, w zależności od warunków w jakich żyją. Zauważono bowiem, że na przestrzeni ostatnich lat ich rozmiar stopniowo maleje i zaczęto się tym faktem niepokoić. Do analizy zebrano pomiary śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek, a w ramach połowu jednej jednostki wybierano losowo od 50 do 100 sztuk trzyletnich śledzi.

Wiersze w zbiorze są uporządkowane chronologicznie, dlatego przyjęto że liczba porządkująca wiersze będzie utożsamiana z czasem.

str(raw_data)
## 'data.frame':    52582 obs. of  16 variables:
##  $ t     : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ length: num  23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
##  $ cfin1 : Factor w/ 40 levels "?","0","0.01",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ cfin2 : Factor w/ 49 levels "?","0","0.01",..: 14 14 14 14 14 14 14 14 14 14 ...
##  $ chel1 : Factor w/ 49 levels "?","0","0.2287",..: 20 20 20 20 20 20 20 20 20 20 ...
##  $ chel2 : Factor w/ 52 levels "?","10.10963",..: NA 23 23 23 23 23 23 23 23 23 ...
##  $ lcop1 : Factor w/ 49 levels "?","0.30741",..: 20 20 20 20 20 20 20 20 20 20 ...
##  $ lcop2 : Factor w/ 52 levels "?","10.72889",..: 23 23 23 23 23 NA 23 23 23 23 ...
##  $ fbar  : num  0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
##  $ recr  : int  482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
##  $ cumf  : num  0.306 0.306 0.306 0.306 0.306 ...
##  $ totaln: num  267381 267381 267381 267381 267381 ...
##  $ sst   : Factor w/ 52 levels "?","12.7690663857",..: 38 38 38 38 38 38 38 38 38 38 ...
##  $ sal   : num  35.5 35.5 35.5 35.5 35.5 ...
##  $ xmonth: int  7 7 7 7 7 7 7 7 7 7 ...
##  $ nao   : num  2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...

Jak widać w rezultacie wywołania powyższej metody, oryginalne, surowe dane zawierają 52582 przykłady opisane przez 16 atrybutów (15 właściwych i 1 porządkujący). Ich dokładniejsza analiza zostanie przedstawiona w dalszej części raportu. W niektorych przypadkach możemy zaobserwować brak pomiaru wartości danego atrybutu oznaczony przez znak zapytania.

Są to zmiene dotyczące dostępności planktonu i temperatury przy powierzchni wody. Dokładne informacje dotyczące brakujących danych z uwzględnieniem kolumn przedstawia poniższa komórka:

## [1] "t: 0 brakujących wartości"
## [1] "length: 0 brakujących wartości"
## [1] "cfin1: 1581 brakujących wartości"
## [1] "cfin2: 1536 brakujących wartości"
## [1] "chel1: 1555 brakujących wartości"
## [1] "chel2: 1556 brakujących wartości"
## [1] "lcop1: 1653 brakujących wartości"
## [1] "lcop2: 1591 brakujących wartości"
## [1] "fbar: 0 brakujących wartości"
## [1] "recr: 0 brakujących wartości"
## [1] "cumf: 0 brakujących wartości"
## [1] "totaln: 0 brakujących wartości"
## [1] "sst: 1584 brakujących wartości"
## [1] "sal: 0 brakujących wartości"
## [1] "xmonth: 0 brakujących wartości"
## [1] "nao: 0 brakujących wartości"

Łącznie daje to nam około 10000 niepełnych przykładów, które zdecydowałyśmy się odfiltrować, aby ułatwić dalszą analizę. Z powodu braku niektórych wartości i oznaczenia ich przez “?” typ niektórych zmiennych został domyślnie ustawiony na Factor (typ wyliczeniowy, może być interpretowany jako kategoria). Widzimy jednak, że pozostałe wartości są numeryczne (całkowite lub niecałkowite), tak więc w dalszej analizie, po odfiltrowaniu niepełnych przykładów, dla ułatwienia zmienimy typ tych zmiennych na numeryczny. R domyślnie ustala typ zmiennych całkowitych na integer, a niecałkowitych na numeric.

dim(data)
## [1] 42488    16

Po odfiltrowaniu brakujących wartości, zbiór danych składa się z 42488 kompletnych przykładów, każdy opisany przez 16 atrybutów. Poniżej wyświetlonych zostało kilka przykładów ze zbioru.

head(data)
##   t length   cfin1   cfin2   chel1    chel2   lcop1    lcop2  fbar   recr
## 2 1   22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 2   25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 3   25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 4   24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 7 6   24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 8 7   23.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
##        cumf   totaln           sst      sal xmonth nao
## 2 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 3 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 4 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 5 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 7 0.3059879 267380.8 14.3069330186 35.51234      7 2.8
## 8 0.3059879 267380.8 14.3069330186 35.51234      7 2.8

Analiza wartości atrybutów

Poniżej przedstawiono analizę wartości poszczególnych atrybutów, które opisują dany rekord. Łącznie danych jest 15 atrybutów. Wśród nich występują zmienne takie jak dostępności różnych gatunków planktonu, informacje dotyczące połowów, właściwości wody etc.

Lista atrybutów jest następująca:

  • length: długość złowionego śledzia [cm];
  • cfin1: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1];
  • cfin2: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2];
  • chel1: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1];
  • chel2: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2];
  • lcop1: dostępność planktonu [zagęszczenie widłonogów gat. 1];
  • lcop2: dostępność planktonu [zagęszczenie widłonogów gat. 2];
  • fbar: natężenie połowów w regionie [ułamek pozostawionego narybku];
  • recr: roczny narybek [liczba śledzi];
  • cumf: łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku];
  • totaln: łączna liczba ryb złowionych w ramach połowu [liczba śledzi];
  • sst: temperatura przy powierzchni wody [°C];
  • sal: poziom zasolenia wody [Knudsen ppt];
  • xmonth: miesiąc połowu [numer miesiąca];
  • nao: oscylacja północnoatlantycka [mb].

Długość złowionego śledzia (length) 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    19.0    24.0    25.5    25.3    26.5    32.5
## [1] "Odchylenie standardowe: 1.649605"

Histogram dystrybucji zmiennej length, jak i postawowe statystyki które ją opisują świadczą o tym że ta zmienna ma rozkład normalny. Średnia oraz media są sobie prawie równe. Sama rozpiętość zmiennej jest stosunkowo mała, najmniejsza zaobserwowana długość to 19 centymetrów, największa - 32.5 centymetra. Należało się tego spodziewać, ponieważ zmienna ta opisuje długość śledzia, a ponadto wiadomo o nich że wszystkie miały około 3 lat. Tę zmienną charakteryzuje także niskie odchylenie standardowe.

Dostępność planktonu - zagęszczenie Calanus finmarchicus 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    2.00   15.00   16.39   28.00   40.00 
## [1] "Odchylenie standardowe: 12.715934"
## [1] "Procent dla którego zmienna osiągała wartość równą minimum: 0.279750"

Na histogramie dystrybucji zmiennej cfin1 widoczna jest zdecydująca przewaga bardzo niskich wartości, bliskich minimum. W prawie 30% przypadkach wartość tej zmiennej wynosiła 2.0. Maksymalne zaobserwowane zagęszczenie wyniosło 40.0. Odchylenie standardowe jest wysokie i rozkład zmiennej nie przypomina żadnego charakterystycznego rozkładu.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   14.00   24.00   23.64   34.00   49.00 
## [1] "Odchylenie standardowe: 13.342136"

W przypadku zmiennej cfin2 także wiele razy zaobserwowano wartości bliskie minimum, jednak nie tak dużo jak w przypadku cfin1. Średnia jak i mediana są w tym przypadku wyższe. Rozpiętość zmiennej jest zbliżona: wartość minimalna wynosi 2.0 a maksymalna 49.0. Zauważalne są trzy słupki histogramu znacznie wyższe od pozostałych. Zmienna ma niestandardowy rozkład i wysokie odchylenie standardowe.

Dostępność planktonu - zagęszczenie Calanus helgolandicus 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   13.00   20.00   23.45   36.00   49.00 
## [1] "Odchylenie standardowe: 13.894213"

Histogram dystrybucji zmiennej chel1 ukazuje, że ta zmienna także ma niestandardowy rozkład. Rozkład ten kumuluje się wokół wartości 13.0, co widać na histogramie jako wysoki słupek. Podstawowe statystyki (w szczególności minimum i maksimum) są zbliżone do tych opisujących zmienną cfin2.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   15.00   27.00   28.09   41.00   52.00 
## [1] "Odchylenie standardowe: 15.388129"

Podobnie dla zmiennej chel2, rozkład jest nieregularny i odchylenie standardowe jest wysokie. Ta zmienna dostępności planktonu ma stosunkowo wysoką średnią. Widoczne są dwa wysokie słupki, tym razem część wartości kumuluje się około liczby 24 (przedział 23-24) i 50 (przedział 49-50).

Dostępność planktonu - zagęszczenie widłonogów 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     2.0    13.0    23.0    23.8    35.0    49.0 
## [1] "Odchylenie standardowe: 13.064853"

Charakterystyka rozkładu zmiennej lcop1 jest podobna jak w przypadku innych zmiennych dotyczących dostępności plakntonu. Zmienna przyjmuje różnorodne wartości; najwięcej wartości przyjmuje w przedziale 22-23, co można zaobserwować w postaci wysokiego słupka na histogramie.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   13.00   27.00   27.83   43.00   52.00 
## [1] "Odchylenie standardowe: 15.769491"

Także w przypadku zmiennej lcop2, jej rozkład jest nieregularny i jej wartości kumulują się w wybranych przedziałach (20-21 oraz 51-52). Odchylenie standardowe jest ponownie wysokie. Zakres wartości tej zmiennej jest podobny jak w przypadku poprzednich zmiennych opisujących dostępność planktonu.

Z uwagi na to, że sześć zmiennych opisujących dostępność planktonu ma podobne charakterystyki (co wynika z tego, co reprezentują), na powyższym wykresie zbiorczo przedstawiono odpowiadające im wykresy pudełkowe. Na wykresie nie są widoczne żadne wartości skrajne, tzw. outliery. Zmienne charakteryzują się podobną rozpiętością, a ich średnie nieznacznie się różnią. Można podejrzewać, że plankton związany z zmienną cfin1 jest zazwyczaj najmniej liczny, a dominującymi gatunkami są chel2 i lcop2.

Natężenie połowów w regionie (fbar) [ułamek pozostawionego narybku] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0680  0.2270  0.3320  0.3306  0.4650  0.8490
## [1] "Odchylenie standardowe: 0.160283"

Średnie natężenie połowów w regionie wynosi 0.3306, co oznacza że średnio około 33% narybku nie zostaje złowione. Z histogramu wynika, że zdarzyły się także (rzadkie) przypadki w których ten współczynnik był bliski maksimum (około 85%). Rozkład zmiennej jest nieregularny, jednak można zaobserwować tendencję to zmniejszania się liczności wystąpień zmiennej im większy ułamek pozostawionego narybku. Większość przypadków nie przekracza 50%.

Roczny narybek (recr) [liczba śledzi] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  140515  360061  421391  519877  724151 1565890
## [1] "Odchylenie standardowe: 270638.842995"

Zmienna recr przyjmuje bardzo wysokie wartości, są to liczby rzędu kilkuset tysięcy. Zaobserwowane maksimum wynosi ponad półtora miliona. Również w tym przypadku odchylenie standardowe jest wysokie. Przeważają wartości rzędu 400 000 śledzi. Zmienna przyjmuje bardzo rożne wartości co jest stosunkowo zaskakujące, biorąc pod uwagę co reprezentuje ta zmienna. Z histogramu wynika, że zdarzały się lata które skutowały w bardzo dużej liczbie narybku, nawet przekraczające średnią dwukrotnie. Można wyróżnić trzy grupy słupków na histogramie, oddzielone przerwami, tj. brakiem wystąpień. Pierwsza z nich jest najbardziej liczna.

Łączne roczne natężenie połowów w regionie (cumf) [ułamek pozostawionego narybku] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.06833 0.14809 0.23191 0.22987 0.29803 0.39801
## [1] "Odchylenie standardowe: 0.092352"

Zmienna cumf jest powiązana ze zmienną fbar - natężenie połowów w regionie, ponieważ dotyczy tej samej obserwacji (ułamku pozostawionego narybku), tylko kumuluje wszystkie wartości zaobserwowane w danym roku. Średnie łączne natężenie połowów wynosi około 23% pozostawionego narybku. Jest to mniej niż wynika ze zmiennej fbar. Można zaobserwować że mimo wystąpienia dla zmiennej fbar wartości skrajnie wysokich, w tym przypadku zakres wartości zmiennej jest już mniejszy. Maksymalne łączne rozczne natężenie połowów w regionie wynosiło niecałe 40%.

Łączna liczba ryb złowionych w ramach połowu (totaln) [liczba śledzi] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  144137  306068  539558  515082  730351 1015595
## [1] "Odchylenie standardowe: 221388.914597"

Podobnie jak zmienna recr, także wyrażana w liczbie śledzi, zmienna totaln osiąga wartości rzędu kilkuset tysięcy, z maksymalną wartością przekraczającą milion. Najwięcej razy ta liczba była rzędu 750 000, co odzwierciedla wysoki słupek na wyrkesie histogramu. Zdarzyły się przypadki połowów w których liczba śledzi była dwa razy większa niż ich średnia liczba.

Temperatura przy powierzchni wody (sst) [°C] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   14.00   25.00   25.45   36.00   52.00
## [1] "Odchylenie standardowe: 13.707105"

Rozkład wartości tej zmiennej jest stosunkowo zaskakujący, ponieważ wynika z niego że często obserwowane były zarówno bardzo niskie jak i bardzo wysokie temperatury wody. Szczególnie maksimum jest zastanawiające - największa wartość tej zmiennej to 52 stopnie celsjusza. Należy brać pod uwagę fakt że pomiar następował przy powierzchni wody, w związku z czym wpływ na taki pomiar z pewnością miało na przykład nasłonecznienie. Średnia temperatura przy powierzchni wody wynosi 25 stopni. Stosunkowo rzadko odnotowano temperatury bliskie zera. Dwa przedziały o największej liczności, jak wynika z histogramu, to 30-35 oraz 10-15.

Poziom zasolenia wody (sal) [Knudsen ppt] 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   35.40   35.51   35.51   35.51   35.52   35.61

Zmienna sal przyjmuje wartości z bardzo małego przedziału, w porównaniu do innych atrybutów tego zbioru danych. Wartości minimalne i maksymalne niewiele różnią się od wartości średniej. Z hsitogramu można wywnioskować, że większość obserwacji przyjęło wartości bliskie średniej (lub mediany, są one sobie równe). Jednak ponadto, część wartości kupuluje się albo przy minimum, albo przy maksimum. Najmniej wystąpień jest wartości pośrednich.

Miesiąc połowu (xmonth) [numer miesiąca] 

## Warning: Ignoring unknown parameters: binwidth, bins, pad

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   5.000   8.000   7.252   9.000  12.000

Najwięcej obserwacji odnotowano dla miesiąca sierpień, w następnej kolejności: październik i lipiec. Najmniej obserwacji jest danych dla stycznia. Warto zwrócić uwagę, że rozkład miesięcy jest nierównomierny. Biorąc dalej pod uwagę, że każdy miesiąc charakteryzują pewne warunki pogodowe, należy uwzględniać taki fakt na przykład przy interpretowaniu średniej temperatury przy powierzchni wody. Wartości z miesiąca sierpień mogłyby zawyżać średnią w związku z tym, że danych jest najwięcej rekordów reprezentujących ten miesiąc. W ogólności z histogramu wynika, że dominują miesiące letnie oraz poprzedzające i następujące ten okres.

Oscylacja północnoatlantycka (nao) [mb] 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -4.89000 -1.90000  0.20000 -0.09642  1.63000  5.08000
## [1] "Odchylenie standardowe: 2.248234"

Zmienna nao jako jedyna przyjmuje wartości ujemne i dodatnie. Jest to zgodne z zakresem wartości jakie przyjmuje tzw. indeks NAO, który może reprezentować wystąpienie dodatniej lub ujemnej fazy NAO. Oscylacja północnoatlantycka jest zjawiskiem meteorologicznym wpływającym na klimat, które manifestuje się poprzez fluktuacje ciśnienia atmosferycznego, temperatury powietrza, prędkości wiatru, ilości opadów i innych parametrów [wikipedia(https://pl.wikipedia.org/wiki/Oscylacja_p%C3%B3%C5%82nocnoatlantycka)]. Na histogramie widać, że rozkład tej zmiennej byłby zbliżony do normalnego, gdyby nie niska liczność wartości w przedziale -2 do -1. Zarówno wartości skrajnie niskie jak i wysokie występują stosunkowo rzadko. Średnia wartość jest bliska zeru. Ta zmienna także ma dość wysokie odchylenie standardowe.

Korelacje między zmiennymi 

corrplot(cor(numeric_data))

Obliczenie wartości współczynnika korelacji pomiędzy zmiennymi, a zwłaszcza korelacji zmiennej length z innymi może zdecydowanie pomóc w odnalezieniu ukrytych zależności prowadzącymi na przykład do odkrycia jakie czynniki mają wpływ na zmniejszenie długości śledzi w ostatnich latach. Współczynnik korelacji przyjmuje wartości od −1 (zupełna korelacja ujemna), przez 0 (brak korelacji) do +1 (zupełna korelacja dodatnia).

Na powyższej macierzy korelacji możemy zaobserwować następujące zależności:

  • silna korelacja dodatnia (około 0.9) pomiędzy zmiennymi fbar oraz cumf. Obywie te zmienne oznaczają ułamek pozostawionego narybku, cumf jest natomiast jego skumulowaną roczną wersją, a więc logicznie, wraz ze wzrostem fbar rośnie cumf tak więc jest to korelacja dodatnia
  • silna korelacja ujemna (około -0.8) pomiędzy zmiennymi totaln i cumf, oraz nieco słabsza (około -0.6) totaln i fbar. totaln oznacza łączną liczbę śledzi złowionych w ramach połowu. Jest to oczywiście wartość odwrotnie proporcjonalna do ułamku pozostawionego narybku w danym czasie.
  • korelacja dodatnia chel1 i lcop1 oraz chel2 i lcop2, które oznaczają zagęszczenie różnych gatunków planktonu. Być może sugeruje to fakt iż wymienione powyżej gatunki kooperują ze sobą i dlatego żyją zagęszczenie działa w sposób pozytywny na zagęszczenie drugich.
  • możemy dostrzec również słabą ujemną zależność długości śledzia od czasu, co potwierdza fakt, iż z czasem ich rozmiar maleje. Widzimy także nieco mocniejszą ujemną zależnośc czasu oraz liczby złowionych śledzi.
  • możemy również dostrzec dość silną ujemną zależność rozmiaru śledzi oraz sst oznaczające temperaturę przy powierzchni wody. Oznacza to że wzraz ze wzrostem temperatury rozmiar śledzi się zmniejsza.
  • natomiast zwiększenie temperatury na powierzchni wody wykazuje pozytywną korelację z nao - oscylacją północno atlantycką. W związku z tym, być może to ona właśnie przyczynia się do jej ocieplenia.

Poniżej dla zobrazowania zależności rozmiaru śledzi od innych czynników przedstawione zostały dwa przykładowe wykresy:

p <- ggplot(data = numeric_data[sample(nrow(numeric_data), 1000), c(2, 13)], aes(x = sst, y = length)) +
  geom_point() +
  geom_smooth()
p
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

  • wykres zależności długości śledzia od temperatury przy powierzchni wody, potwierdza obserwację z macierzy korelacji - widać trend, że wzraz ze wzrostem temperatury rozmiar śledzi maleje.
p <- ggplot(data = numeric_data[sample(nrow(numeric_data), 1000), c(2, 10)], aes(x = recr, y = length)) +
  geom_point() +
  geom_smooth()
p
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

  • wykres zależności długości śledzia od liczby śledzi, również potwierdza obserwację z macierzy korelacji - nie widać żadnej konkretnej zależności.

Zmiana rozmiaru śledzi w czasie 

length_in_time <- raw_data[sample(nrow(raw_data), 1000), 1:2]
length_in_time <- data.frame(sapply(length_in_time, as.numeric))

p <- ggplot(data = length_in_time, aes(x = t, y = length)) +
  geom_point() +
  geom_smooth()

p <- ggplotly(p)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p

Powyższy interaktywny wykres przedstawia zależność rozmiaru śledzi od czasu. Z względów zasobowych pobrana została próbka 1000 przykładów spośród wszystkich 52582 (aby uzyskać lepsze uogólnienie, było to możliwe, ponieważ każdy rekord zawierał pomiar długości śledzia). Aby wyświetlić dokładną wartość interesującego użytkownika punktu, może on najechać myszką na dany obszar i odczytać wartość.

Smoothed conditional means, czyli wygładzona średnia długości śledzia w czasie została dodana do wykresu, aby łatwiej było zauważyć wzorce występujące w danych. Do wygładzania została użyta domyślna metoda gam oraz formuła y ~ s(x, bs = “cs”). Poświata która znajduje się pod wygładzoną średnią oznacza niski oraz wysoki punktowy przedział ufności wokół średniej.

Dzięki linii trendu na wykresie można z łatwością zaobserować, że do około 15000 przykładu rozmiar śledzia rósł, a później zaczął stopniowo spadać, aż osiągnął średnie minimum w ostatnim przykładzie.

Regresor przewidujący rozmiar śledzia

Przygotowanie zbioru trenującego i testowego

Aby poprawnie dobrać parametry modelu oraz oszacować jego skuteczność, zbiór danych został podzielony na dane uczące, walidujące i testowe w stosunku 0.56 do 0.24 do 0.20. W rezultacie otrzymane zostały następujące rozmiary zbiorów:

## [1] "Rozmiary zbiorów:"
## [1] "Trenujący: 23793 | Walidacyjny: 10197 | Testowy: 8498"

Regresja liniowa

Podstawową metodą wykorzystywaną przy zadaniu regresji jest regresja liniowa. Można do tego wykorzystać implementację takiego modelu o nazwie lm. Bez wątpienia zaletą regresji liniowej jest prostota, łatwość interpretacji wyników i małe zapotrzebowanie na moc obliczeniową.

Poniżej zaprezentowano parametry modelu otrzymane po wytrenowaniu go na zbiorze trenującym.

## 
## Call:
## lm(formula = length ~ ., data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1392 -0.8885  0.0499  0.9098  6.1734 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.423e+02  9.347e+00  15.222  < 2e-16 ***
## cfin1        3.007e-02  9.306e-04  32.313  < 2e-16 ***
## cfin2       -5.657e-03  9.588e-04  -5.900 3.68e-09 ***
## chel1        1.457e-02  1.126e-03  12.941  < 2e-16 ***
## chel2        1.546e-02  9.791e-04  15.790  < 2e-16 ***
## lcop1       -3.968e-03  1.222e-03  -3.246  0.00117 ** 
## lcop2       -1.095e-02  1.019e-03 -10.738  < 2e-16 ***
## fbar         5.405e+00  1.260e-01  42.885  < 2e-16 ***
## recr        -7.646e-07  4.158e-08 -18.388  < 2e-16 ***
## cumf        -9.360e+00  2.666e-01 -35.108  < 2e-16 ***
## totaln      -3.583e-07  7.200e-08  -4.976 6.54e-07 ***
## sst         -4.526e-02  9.169e-04 -49.356  < 2e-16 ***
## sal         -3.258e+00  2.635e-01 -12.367  < 2e-16 ***
## xmonth       1.262e-02  3.166e-03   3.985 6.76e-05 ***
## nao          6.799e-03  6.316e-03   1.077  0.28169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.334 on 23778 degrees of freedom
## Multiple R-squared:  0.3413, Adjusted R-squared:  0.341 
## F-statistic: 880.2 on 14 and 23778 DF,  p-value: < 2.2e-16

Z wywołania metody summary można odczytać, które zmienne okazały się znaczące. Na etapie poszukiwania najlepszego bardziej istotne są jednak wartości Residual standard error oraz R-squared.

Prosty model regresji liniowej wykorzystujący wszystkie atrybuty jakimi można byłoby opisać zmienną length myli się średnio o 1.34 centymetry w stosunku do prawdziwej długości śledzia.

Regresja liniowa z reularyzacją (Ridge, Lasso i Elastic-Net)

Regresję liniową można ulepszyć stosując regularyzację. W języku R dostępna jest paczka o nazwie glmnet, która pozwala na przetestowanie w prosty sposób różne rodzaje regresji z regularyzacją. W raporcie zaprezentowane zostaną 3 z nich:

  • Ridge, stosowana aby zmniejszyć problem współliniowości w regresji liniowej, zapewniając ulepszoną wydajność w problemie estymacji parametrów. Model jest mniej podatny na szum w danych, za to zmniejszona jest jego zdolność dopasowania do danych. W funkcji straty dodatkowo występuje kara L2.

  • Lasso, używana do polepszenia dokładności predykcji poprzez wybór zmiennych dla niej istotnych (również kosztem mniejszego dopasowania do danych). W funkcji straty dodatkowo występuje kara L1.

  • Elastic-Net, w której występuje liniowa kombinacja kar L1 oraz L2 z wymienionych wcześniej metod.

Poprzez użycie funkcji cv.glmnet z parametrami type.measure = "mse" oraz family = "gaussian" będziemy trenować model regresji liniowej. Jednak oprócz standardowego minimalizowania sumy kwadratów odchyłek dodana jest także kara w postaci: \[\lambda *[\alpha*(|var_1| + ... + |var_n|) + (1-\alpha)*(var_1^2 + ... + var_n^2)]\] Wartość zmiennej \(\lambda\) dobierana jest metodą kroswalidacji. Pierwszy człon powyższej formuły odpowiada za karę L1 pomnożoną przez współczynnik \(\alpha\), drugi za karę L2 pomnożoną przez współczynnik \(1-\alpha\). Dzięki temu dobierając różne wartości parametru \(\alpha\) jesteśmy w stanie przetestować zarówno regularyzację Ridge, jak i Lasso oraz Elastic-Net. Dla \(\alpha=0\) otrzymujemy Ridge, dla \(\alpha=1\) Lasso, a dla pozostałych wartości \(\alpha\) Elastic-Net.

##    alpha     rmse       rsq fit.name
## 1    0.0 1.362720 0.3298651   alpha0
## 2    0.1 1.356963 0.3355157 alpha0.1
## 3    0.2 1.358662 0.3338504 alpha0.2
## 4    0.3 1.359021 0.3334983 alpha0.3
## 5    0.4 1.358608 0.3339031 alpha0.4
## 6    0.5 1.357837 0.3346588 alpha0.5
## 7    0.6 1.357347 0.3351390 alpha0.6
## 8    0.7 1.359317 0.3332078 alpha0.7
## 9    0.8 1.358077 0.3344237 alpha0.8
## 10   0.9 1.357830 0.3346661 alpha0.9
## 11   1.0 1.357638 0.3348537   alpha1
p <- ggplot() +
  geom_point(data = results, aes(x = alpha, y = rmse))
p

p <- ggplot() +
  geom_point(data = results, aes(x = alpha, y = rsq))
p

Na powyższych wykresach zaprezentowane zostały wartości \(RMSE\) oraz \(R-squared\) otrzymane dla zbioru walidacyjnego, po wytrenowaniu modelu o różnych wartościach parametru \(\alpha\) na zbiorze trenującym. Jak można zauważyć, najmniejszą wartość \(RMSE\) oraz najlepsze dopasowanie \(R-squared\) uzyskane zostało dla współczynnika \(\alpha=0.4\).

Regresja nieliniowa (wielomianowa)

Często okazuje się, że zależność między zmiennymi niestety nie jest liniowa, a wtedy zwykła regresja liniowa może okazać się niewystarczająca. W takim przypadku można jednak spróbować zastosować regresję wielomianową.

W R jest możliwość zaprojektowania takiej metody poprzez dodanie zmiennych będących istniejącymi zmiennymi podniesionymi do potęgi. W naszym przypadku, po przetestowaniu różnych wariantów, najlepszy wynik uzyskałyśmy dodając zmienne podniesione do 2 oraz do 3 potęgi.

Poniżej zaprezentowano parametry modelu otrzymane po wytrenowaniu go na zbiorze trenującym.

lm_mod <- lm(length ~ poly(cfin1, 3) + poly(cfin2, 3) + poly(chel1, 3) + poly(chel2, 3)  + poly(lcop1, 3)  + poly(lcop2, 3) + poly(fbar, 3) + poly(recr, 3) + poly(cumf, 3) + poly(totaln, 3) + poly(sst, 3) + poly(sal, 3) + poly(xmonth, 3) + poly(nao, 3) , data = train)

summary(lm_mod)
## 
## Call:
## lm(formula = length ~ poly(cfin1, 3) + poly(cfin2, 3) + poly(chel1, 
##     3) + poly(chel2, 3) + poly(lcop1, 3) + poly(lcop2, 3) + poly(fbar, 
##     3) + poly(recr, 3) + poly(cumf, 3) + poly(totaln, 3) + poly(sst, 
##     3) + poly(sal, 3) + poly(xmonth, 3) + poly(nao, 3), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8653 -0.8347  0.0143  0.8396  5.8613 
## 
## Coefficients:
##                    Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)       2.530e+01  8.079e-03 3130.980  < 2e-16 ***
## poly(cfin1, 3)1   1.866e+01  2.748e+00    6.792 1.13e-11 ***
## poly(cfin1, 3)2  -3.163e+01  2.220e+00  -14.248  < 2e-16 ***
## poly(cfin1, 3)3   3.054e+01  1.745e+00   17.502  < 2e-16 ***
## poly(cfin2, 3)1   2.976e+01  3.670e+00    8.109 5.33e-16 ***
## poly(cfin2, 3)2  -2.166e+01  2.229e+00   -9.721  < 2e-16 ***
## poly(cfin2, 3)3   1.777e+00  2.099e+00    0.846  0.39731    
## poly(chel1, 3)1   2.036e+01  4.114e+00    4.949 7.51e-07 ***
## poly(chel1, 3)2   5.253e+01  3.671e+00   14.311  < 2e-16 ***
## poly(chel1, 3)3   1.511e+01  3.287e+00    4.598 4.30e-06 ***
## poly(chel2, 3)1  -2.439e+01  4.666e+00   -5.227 1.74e-07 ***
## poly(chel2, 3)2  -1.287e+02  5.292e+00  -24.318  < 2e-16 ***
## poly(chel2, 3)3   1.964e+00  3.019e+00    0.650  0.51548    
## poly(lcop1, 3)1   1.894e+01  4.094e+00    4.627 3.74e-06 ***
## poly(lcop1, 3)2   1.159e+01  3.786e+00    3.062  0.00220 ** 
## poly(lcop1, 3)3  -2.095e+01  2.992e+00   -7.002 2.59e-12 ***
## poly(lcop2, 3)1   3.763e+01  5.550e+00    6.781 1.22e-11 ***
## poly(lcop2, 3)2   9.425e+01  5.026e+00   18.754  < 2e-16 ***
## poly(lcop2, 3)3  -2.469e+01  3.880e+00   -6.365 1.99e-10 ***
## poly(fbar, 3)1    2.515e+02  8.281e+00   30.369  < 2e-16 ***
## poly(fbar, 3)2   -8.618e+01  4.999e+00  -17.240  < 2e-16 ***
## poly(fbar, 3)3    1.620e+01  3.560e+00    4.550 5.39e-06 ***
## poly(recr, 3)1   -9.732e+00  3.262e+00   -2.983  0.00286 ** 
## poly(recr, 3)2    1.758e+01  2.209e+00    7.958 1.82e-15 ***
## poly(recr, 3)3   -3.472e+01  2.441e+00  -14.224  < 2e-16 ***
## poly(cumf, 3)1   -3.107e+02  1.031e+01  -30.139  < 2e-16 ***
## poly(cumf, 3)2    7.751e+01  4.733e+00   16.377  < 2e-16 ***
## poly(cumf, 3)3   -3.875e+01  2.653e+00  -14.604  < 2e-16 ***
## poly(totaln, 3)1 -6.387e+01  4.934e+00  -12.945  < 2e-16 ***
## poly(totaln, 3)2 -3.201e+00  4.241e+00   -0.755  0.45044    
## poly(totaln, 3)3  3.051e+01  3.696e+00    8.255  < 2e-16 ***
## poly(sst, 3)1    -7.321e+01  3.322e+00  -22.038  < 2e-16 ***
## poly(sst, 3)2    -7.131e+01  3.000e+00  -23.770  < 2e-16 ***
## poly(sst, 3)3     1.484e+01  2.501e+00    5.933 3.02e-09 ***
## poly(sal, 3)1    -7.356e+01  2.603e+00  -28.255  < 2e-16 ***
## poly(sal, 3)2     7.950e+01  3.709e+00   21.432  < 2e-16 ***
## poly(sal, 3)3    -3.360e+01  3.134e+00  -10.724  < 2e-16 ***
## poly(xmonth, 3)1  3.085e+00  1.274e+00    2.422  0.01544 *  
## poly(xmonth, 3)2 -2.037e+01  1.266e+00  -16.092  < 2e-16 ***
## poly(xmonth, 3)3 -6.189e+00  1.272e+00   -4.867 1.14e-06 ***
## poly(nao, 3)1    -2.978e+01  3.303e+00   -9.016  < 2e-16 ***
## poly(nao, 3)2    -2.439e+01  2.796e+00   -8.722  < 2e-16 ***
## poly(nao, 3)3     4.523e+01  2.339e+00   19.332  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.246 on 23750 degrees of freedom
## Multiple R-squared:  0.4261, Adjusted R-squared:  0.4251 
## F-statistic: 419.8 on 42 and 23750 DF,  p-value: < 2.2e-16

Z wywołania funkcji summary możemy wywnioskować, że niektóre atrybuty okazały się nieistotne. Są to na przykład: zmienna cfin2 podniesiona do trzeciej potęgi oraz zmienna chel2 również podniesiona do potęgi trzeciej. W dalszej analizie można by usunąć te zmienne z modelu w celu jego ulepszenia.

Poniżej przedstawione zostały wartości parametrów \(RMSE\) oraz \(R-squared\) otrzymane dla zbioru walidacyjnego, po wytrenowaniu modelu na zbiorze trenującym.

pred <- predict(lm_mod, validation)
printf("RMSE: %f, R2: %f",  calculate_rmse(pred, validation$length, length(pred)), calculate_rsq(pred, validation$length))
## [1] "RMSE: 1.256971, R2: 0.429836"

Jak widać, model regresji wielomianowej wykorzystujące wszystkie atrybuty podniesione do 1, 2 oraz 3 potęgi myli się średnio o 1.26 centymetra w stosunku do prawdziwej długości śledzia, co daje nam lepszy wynik niż w przypadku regresji liniowej. Zgodnie z oczekiwaniami, model cechuje się lepszym dopasowaniem do danych, ze względu na fakt większego stopnia skoplikowania zależności między zmiennymi w rozpatrywanym przez nas problemie. W związku z tym faktem postanowiłyśmy również przetestować algorytm bardziej dla niego odpowiedni.

Gradient Boosting Machine

Kolejnym algorytmem który można wykorzystać do zadania regresji jest Gradient Boosting Machine. Jest to algorytm z rodziny tak zwanych ensemble methods, metod wykorzystujących wiele algorytmów uczenia celem uzyskania trafniejszych wyników.

W rozważanym problemie skorzystałyśmy z kombinacji drzew decyzyjnych, co jest często spotykanym rozwiązaniem w kontekście Gradient Boosting Machine. Takie rozwiązanie jest zaimplementowane w bibliotece gbm.

Liczba drzew generowanych w modelu wynosi 2000. Każde drzewo opisuje parametr interaction.depth, który ma wartość 4 i określa że każde drzewo jest małym drzewem z 4 podziałami. Parametr shrinkage jest swego rodzaju prędkością uczenia (learning rate). Model wykorzystuje wszystkie atrybuty, jakie są dostępne.

NUM_TREES <- 2000
model <- gbm(length ~ . ,data = train,distribution = "gaussian",n.trees = NUM_TREES, shrinkage = 0.1, interaction.depth = 4)
model
## gbm(formula = length ~ ., distribution = "gaussian", data = train, 
##     n.trees = NUM_TREES, interaction.depth = 4, shrinkage = 0.1)
## A gradient boosted model with gaussian loss function.
## 2000 iterations were performed.
## There were 14 predictors of which 14 had non-zero influence.

Po wytrenowaniu modelu możliwa jest analiza jego skuteczności na zbiorze walidacyjnym. W tym szczególnym przypadku można także zaobserwować jak zmniejsza się błąd kwadratowy w miarę jak dodawane są nowe drzewa (zgodnie z koncepcją boostingu, nowe drzewa “uczą się na błędach poprzednich”).

n_trees_seq = seq(from=100, to=NUM_TREES, by=100)
pred <- predict(model, validation, n.trees=n_trees_seq)

mses <- with(validation, apply( (pred - length) ^ 2, 2, mean))
rmses <- sqrt(mses)

rss <- with(validation, apply( (pred - length) ^ 2, 2, sum))
tss <- sum((y_val_var - mean(y_val_var)) ^ 2)
rsq <- 1 - rss/tss

Poniżej zaprezentowano początkowe wartości błędu \(RMSE\) na zbiorze walidacyjnym:

##      100      200      300      400      500      600 
## 1.182680 1.168758 1.162520 1.159144 1.156600 1.154503

A także końcowe, aż do osiągnięcia liczby drzew równej 2000:

##     1500     1600     1700     1800     1900     2000 
## 1.150461 1.150446 1.149929 1.150144 1.150111 1.150460

Dodawanie kolejnych drzew pozwala zmniejszyć błąd predykcji długości śledzia aż do 1.1327.

Poniższa komórka wypisuje ponadto osiągnięte wartości \(R^2\) dla końcowych iteracji dla których liczba drzew była największa.

##      1500      1600      1700      1800      1900      2000 
## 0.5223690 0.5223809 0.5228104 0.5226320 0.5226593 0.5223693

Analiza ważności atrybutów najlepszego modelu (Gradient Boosting Machine) <a name="gbm_analysis>

Modelem dla którego uzyskano najlepsze wyniki na zbiorze walidacyjnym okazał się model Gradient Boosting Machine z drzewami decyzyjnymi. W związku z tym, ten model został wykorzystany do predykcji długości śledzia w oparciu o dane atrybuty na zbiorze testowym.

Poniżej zaprezentowano kod wywołujący predykcję na zbiorze testowym i wyliczający miary oceny tego modelu na tym zbiorze.

y_test_var <- test$length

n_trees_seq = seq(from=100, to=NUM_TREES, by=100)
pred <- predict(model, test, n.trees=n_trees_seq)

mses <- with(test, apply( (pred - length) ^ 2, 2, mean))
rmses <- sqrt(mses)

rss <- with(test, apply( (pred - length) ^ 2, 2, sum))
tss <- sum((y_test_var - mean(y_test_var)) ^ 2)
rsq <- 1 - rss/tss

Wartości \(RMSE\) dla ostatnich iteracji - modelu z największą liczbą drzew - zaprezentowano poniżej:

##     1500     1600     1700     1800     1900     2000 
## 1.138825 1.138008 1.138558 1.137338 1.136831 1.136998

Podobnie dla miary \(R^2\):

##      1500      1600      1700      1800      1900      2000 
## 0.5225076 0.5231925 0.5227314 0.5237534 0.5241781 0.5240379

Redukcję błędu \(MSE\) na zbiorze testowym wraz z wzrostem liczby drzew przedstawia poniższy wykres:

plot(n_trees_seq, mses, pch=19, col="blue", xlab="Liczba drzew", ylab="MSE", main="Efektywność boostingu na zbiorze testowym")

Analizę ważności atrybutów wykorzystanych do przewidzenia długości śledzia można przeprowadzić z użyciem funkcji summary.

##           var   rel.inf
## sst       sst 36.973848
## xmonth xmonth 10.242841
## recr     recr  9.755441
## totaln totaln  8.767468
## cfin2   cfin2  5.100975
## lcop1   lcop1  4.768724
## fbar     fbar  4.513241
## cfin1   cfin1  3.693603
## nao       nao  3.626051
## sal       sal  2.822131
## chel1   chel1  2.820460
## chel2   chel2  2.504430
## lcop2   lcop2  2.382701
## cumf     cumf  2.028087

Metodą ustalania ważności atrybutów wykorzystaną w tej funkcji (z takim modelem jak gbm) jest metoda relative.influence, czyli względny wpływ. W ten sposób, najbardziej istotne są zmienne o największym średnim spadku błędu średniokwadratowego (\(MSE\)), obliczanego przy każdym podziale w każdym drzewie.

Z wykresu wynika, że najistotniejszą zmienną jest sst odpowiadająca za temperaturę przy powierzchni wody. Z uprzedniej analizy wiadomo, że długość śledzia i temperatura przy powierzchni wody są negatywnie skorelowane, co oznacza że im wyższa temperatura tym krótsze są śledzie.

Model uznał także zmienną xmonth za stosunkowo istotną mimo że nie była ona skorelowana z długością śledzia w żaden sposób. Równie ważne są zmienne reprezentujące roczny narybek (recr) i łączną liczbę ryb złowionych w ramach połowu (totaln). Zmienne dotyczące planktonu, potencjalnego pożywienia czy w pewien sposób wskaźnika jakości wody, są stosunkowo nieważne, podobnie jak zasolenie wody.

Do analizy ważności atrybutów modelu można wykorzystać także bibliotekę lime. Pozwala ona na przyjrzenie się pojedynczym obserwacjom i uzasadnienie otrzymanej wartości predykcji. W porównaniu do pierwszej wizualizacji, ukazującej średnie wpływy zmiennych na wynik, poniższe wykresy są dopasowane z dokładnością do jednego przykładu.

Powyższe wykresy są zróżnicowane. Powtarza się w nich odwołanie do zmiennej sst. Bardzo silny wpływ ma zmienna sst na długość śledzia kiedy temperatura przy powierzchni wody jest bardzo wysoka - patrz Case 36958, wpływ ten jest wówczas negatywny. Z kolei pozytywny wpływ na długość śledzia zaobserwowano gdy temperatura nie przekraczała 14 stopni Celsjusza - Case 29236. W jednym z pokazanych przypadków największe znaczenie miała zmienna recr, a zaraz po niej istotna ponownie okazała się sst.

Bardziej skoncentrowaną reprezentację można uzyskać za pomocą metody plot_explanations, która pozwoli zwizualizować więcej przypadków jednocześnie.

## Warning: Unknown or uninitialised column: 'label'.

Na tzw. mapie cieplnej (heat map) zaprezentowano 20 losowo wybranych przypadków z zbioru trenującego. Znaczący wpływ temperatury przy powierzchni o wysokiej wartości, przekraczającej 37 stopni, powtarza się w wielu przypadkach. Podobnie jest z temperaturą niską. Umiarkowany negatywny wpływ na długość śledzia ma fbar <= 0.227, co oznaczałoby że niskie natężenie połowów w regionie wiąże się z zmniejszeniem długości śledzia.

Należy przy okazji tej analizy podkreślić, że najlepsze wyniki udało się uzyskać implementując najbardziej skomplikowany z wykorzystywanych modeli. Bez wątpienia ma on większe zdolności do modelowania skomplikowanych zależności niż podstawowe modele regresji. Nie można jednak mówić tutaj o zbyt skomplikowanym modelu (overfitting), ponieważ poradził on sobie równie dobrze na zbiorze testowym który został użyty tylko raz, na koniec analizy.